///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/ham/tei.h
///
/// \brief   header file for two electron integral functions
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef FORMIC_HAM_TEI_HEADER
#define FORMIC_HAM_TEI_HEADER

namespace formic {

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   Returns the number of elements in a compact two electron integral array.
  ///          Such an array takes advantage of the 8-fold symmetry of the two electron integrals.
  /// 
  /// \param[in]    n           number of orbitals in the basis
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  inline int tei_compact_array_length(const int n) {
    const int nn1 = n * ( n + 1 );
    return nn1 * ( nn1 + 2 ) / 8;
  }

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   compute the compound index for accessing a two electron integral
  /// 
  /// \param[in]  p   1st index for electron coordinate 1
  /// \param[in]  q   2nd index for electron coordinate 1
  /// \param[in]  r   1st index for electron coordinate 2
  /// \param[in]  s   2nd index for electron coordinate 2
  ///
  /// \return  the compound index pqrs
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  inline int tei_compound_index(const int p, const int q, const int r, const int s) {

    // define some local functions that help us avoid branching
    struct LocalFuncs {

      static int max(const int a, const int b) {
        return ( a - b > 0 ? a : b );
      }

      static int cmpd(const int p, const int q) {
        const int px = LocalFuncs::max(p,q);
        const int qx = p + q - px;
        return px * ( px + 1 ) / 2 + qx;
      }

    };

    // compute the compound index
    const int pq = LocalFuncs::cmpd(p,q);
    const int rs = LocalFuncs::cmpd(r,s);
    const int pqrs = LocalFuncs::cmpd(pq,rs);

    // return the total compound index
    return pqrs;

  }

  template<class S> void tei_compact_to_dense(const int n, const S * const tei_cmp, S * const tei_dns);
  template<class S> void tei_dense_to_compact(const int n, const S * const tei_dns, S * const tei_cmp);
  template<class S> void tei_rotate_indices_left(const int n, const S * const tei_in, S * const tei_out);
  template<class S> void tei_rotate_indices_right(const int n, const S * const tei_in, S * const tei_out);
  template<class S> void tei_swap_first_two_indices(const int n, const S * const tei_in, S * const tei_out);
  template<class S> void tei_swap_last_two_indices(const int n, const S * const tei_in, S * const tei_out);
  template<class S> void tei_get_ovov_aaaa_slice(const int no,
                                                 const int nv,
                                                 const int * const ao2s,
                                                 const int * const bo2s,
                                                 const int * const au2s,
                                                 const int * const bu2s,
                                                 const S * const tei_dense,
                                                       S * const tei_slice);
  template<class S> void tei_get_ovov_aabb_slice(const int no,
                                                 const int nv,
                                                 const int * const ao2s,
                                                 const int * const bo2s,
                                                 const int * const au2s,
                                                 const int * const bu2s,
                                                 const S * const tei_dense,
                                                       S * const tei_slice);

}

#endif
